Abstract
Part of theR for Artists and Designers course at the School of Foundation Studies, Srishti Manipal Institute of Art, Design, and Technology, Bangalore.
At the end of this Lab session, we should: - know the types and structures of spatial data and be able to work with them - understand the basics of modern spatial packages in R - be able to specify and download spatial data from the web, using R - plot static and interactive maps using ggplot, tmap and leaflet packages - add symbols and markers for places and regions of our own interest in these maps. - see directions for further work (e.g. maps + networks together)
We will take small steps in making maps using just two of the several map making packages in R.
The steps we will use are:
prettymapr or similar..)osmdataosmplot, tmap and also with leaflet.leaflet using a variety of map data providers. Note: tmap can also do interactive maps which we will explore also.Bas. Onwards and Map-wards!!
Let’s get BLR data into R and see if we can plot an area of interest. Then we can order on Swiggy and…never mind.
Where is my home? Specify a “bounding box” first, using a rough longitude latitude info directly, or using a place name to search for the long/lat info:
# BLR Bounding Box
# get_bbox needs lat and lon ranges
bbox <- osmplotr::get_bbox(c(77.56,12.93,77.63,12.96))
bbox_l <- osmdata::getbb("Bangalore, India") # LARGE
bbox_p <- prettymapr::searchbbox("Bangalore") # ALSO LARGE
## Using default API key for pickpoint.io. If batch geocoding, please get your own (free) API key at https://app.pickpoint.io/sign-up
bbox
## min max
## x 77.56 77.63
## y 12.93 12.96
bbox_l
## min max
## x 77.46010 77.78405
## y 12.83401 13.14366
bbox_p # identical with bbox_l
## min max
## x 77.46010 77.78405
## y 12.83401 13.14366
# Get Map data
dat_B <- extract_osm_objects (key = "building", bbox = bbox)
dat_H <- extract_osm_objects (key = 'highway', bbox = bbox)
dat_P <- extract_osm_objects (key = 'park', bbox = bbox)
dat_G <- extract_osm_objects (key = 'landuse', value = 'grass', bbox = bbox)
dat_T <- extract_osm_objects (key = 'natural', value = 'tree', bbox = bbox)
# Useful keys include building, place, amenity, shop, waterway, natural, boundary, and highway.
# How many buildings?
nrow(dat_B)
## [1] 39699
dat_B$geometry
## Geometry set for 39699 features
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 77.55967 ymin: 12.92885 xmax: 77.63028 ymax: 12.96168
## Geodetic CRS: WGS 84
## First 5 geometries:
## POLYGON ((77.58405 12.93005, 77.5845 12.93005, ...
## POLYGON ((77.58647 12.94997, 77.58645 12.94994,...
## POLYGON ((77.61162 12.9368, 77.61159 12.93667, ...
## POLYGON ((77.61255 12.93615, 77.61251 12.93601,...
## POLYGON ((77.61291 12.93646, 77.61276 12.93594,...
class(dat_B$geometry)
## [1] "sfc_POLYGON" "sfc"
So dat_B has 39699 buildings and their geometry is naturally a POLYGON type of geometry column.
Do this check for all the other spatial data.
We could quickly plot this using the package osmplotr. However, in my ( i.e. Arvind’s opinion ) it is not as flexible as other packages. Maybe I need to study it in more detail.
So we will continue with ggplot:
blr_map <- ggplot() +
geom_sf(data = dat_B, colour = "orange") + # POLYGONS
geom_sf(data = dat_H, col = "gray20") + # LINES
geom_sf(data = dat_G, col = "darkseagreen1") +
geom_sf(data = dat_P, col = "darkseagreen") +
geom_sf(data = dat_T, col = "green") # POINTS
blr_map
Note how geom_sf is capable of handling any geometry in the sfc column !! > geom_sf() is an unusual geom because it will draw different geometric objects depending on what simple features are present in the data: you can get points, lines, or polygons.
We can create areas of interest around the map. For a start: we are simply going to zoom in to an area and say that is our area of interest. This area then needs to go Through the Looking Glass and become part of the projection of dat_B!!
We can then try based on your favourite restaurants etc.
my_area <- bbox %>%
# zooming in within our bounding box area
prettymapr::zoombbox(factor = 8, offset = c(0,0))
my_area
## min max
## x 77.59063 77.59937
## y 12.94313 12.94688
bbox
## min max
## x 77.56 77.63
## y 12.93 12.96
OK, my area is smalled and contained inside my bbox. So zoom works. Now to convert my_area into a spatial dataframe using sf and add it to the blr_map plot:
my_area_in_blr <- matrix(
c(
my_area["x", "min"],
my_area["x", "max"],
my_area["x", "max"],
my_area["x", "min"],
my_area["x", "min"],
my_area["y", "min"],
my_area["y", "min"],
my_area["y", "max"],
my_area["y", "max"],
my_area["y", "min"]),
ncol = 2) %>% # Make a matrix
list() %>% # Convert to list since POLYGON needs a list
st_polygon() %>% # Convert to POLYGON
st_sfc(., crs = st_crs(dat_B)) %>% # Convert POLYGON to an `sfc`
# Through the Looking Glass with dat_B
st_as_sf(.) # Convert sfc to an sf spatial dataframe. Phew!!
my_area_in_blr
## Simple feature collection with 1 feature and 0 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 77.59062 ymin: 12.94312 xmax: 77.59937 ymax: 12.94688
## Geodetic CRS: WGS 84
## x
## 1 POLYGON ((77.59063 12.94313...
blr_map <-
blr_map + geom_sf(data = my_area_in_blr, colour = "red")
blr_map
Is it time to order on Swiggy…
Let us adding interesting places to our map: say based on your favourite restaurants etc. We need restaurant data: lat/long + name + maybe type of restaurant. This can be manually created ( like all of OSMdata ) or if it is already there we can download using key-value pairs in our OSM data query.
dat_R <- extract_osm_objects(bbox = bbox, key = "amenity", value = "restaurant", return_type = "point")
## Issuing query to Overpass API ...
## Rate limit: 0
## Query complete!
## converting OSM data to sf format
# How many restaurants have we got?
dat_R %>% nrow()
## [1] 300
names(dat_R)
## [1] "osm_id" "name" "addr.city"
## [4] "addr.country" "addr.floor" "addr.hamlet"
## [7] "addr.housename" "addr.housenumber" "addr.place"
## [10] "addr.postcode" "addr.street" "alt_name"
## [13] "amenity" "capacity" "contact.phone"
## [16] "cuisine" "delivery" "description"
## [19] "designation" "diet.vegan" "diet.vegetarian"
## [22] "email" "food" "internet_access"
## [25] "level" "name.en" "name.kn"
## [28] "note" "opening_hours" "operator"
## [31] "outdoor_seating" "phone" "smoking"
## [34] "source" "speciality" "takeaway"
## [37] "toilets.wheelchair" "website" "wheelchair"
## [40] "wikidata" "wikipedia" "geometry"
# Let's look at the `cuisine` column!
# ( I want pizza...)
dat_R$cuisine
## [1] NA NA
## [3] NA NA
## [5] NA NA
## [7] NA "indian"
## [9] NA NA
## [11] "indian" NA
## [13] "chinese" NA
## [15] NA NA
## [17] NA NA
## [19] "regional" NA
## [21] "Mughlai" NA
## [23] "chinese" "north indian"
## [25] NA "chinese"
## [27] "chinese" NA
## [29] "indian" NA
## [31] "tandoor" NA
## [33] NA NA
## [35] "regional" "Naga"
## [37] "chinese" NA
## [39] NA NA
## [41] "chinese" NA
## [43] NA NA
## [45] NA "indian"
## [47] "spanish" NA
## [49] "indian" "chinese"
## [51] NA NA
## [53] NA "indian"
## [55] NA "indian"
## [57] "indian" "ice_cream"
## [59] "chinese" "indian"
## [61] NA "ice_cream"
## [63] "indian" NA
## [65] "indian" "indian"
## [67] "international" NA
## [69] "international" NA
## [71] "ice_cream" "seafood"
## [73] NA NA
## [75] "indian" "international"
## [77] "chinese" "indian"
## [79] "indian" "chinese"
## [81] "indian" "international"
## [83] "ice_cream" "ice_cream"
## [85] NA "italian"
## [87] NA NA
## [89] "indian" "chinese"
## [91] "ice_cream" "international"
## [93] "indian" NA
## [95] NA "coffee_shop"
## [97] "international" NA
## [99] NA NA
## [101] "indian" "international"
## [103] "international" "international"
## [105] "international" "indian"
## [107] NA NA
## [109] NA NA
## [111] NA NA
## [113] "indian" "international"
## [115] NA NA
## [117] NA NA
## [119] "chinese" NA
## [121] "indian" NA
## [123] NA "indian"
## [125] "chinese" "indian"
## [127] "indian" "indian"
## [129] NA "international"
## [131] "indian" NA
## [133] "italian" NA
## [135] "indian" "indian"
## [137] "regional" "indian"
## [139] "international" NA
## [141] "regional" NA
## [143] NA NA
## [145] "chinese" "pizza"
## [147] NA NA
## [149] NA NA
## [151] NA NA
## [153] NA NA
## [155] NA NA
## [157] NA NA
## [159] NA NA
## [161] NA NA
## [163] "pizza" "Multi-cuisne"
## [165] "indian" "indian"
## [167] "indian" NA
## [169] "andhra" NA
## [171] NA NA
## [173] NA NA
## [175] "regional" NA
## [177] "regional" "regional"
## [179] "regional" "regional"
## [181] "regional" "regional"
## [183] "regional" "regional"
## [185] "regional" "regional"
## [187] "regional" "regional"
## [189] "regional" "regional"
## [191] "regional" "regional"
## [193] "regional" "regional"
## [195] "italian" "regional"
## [197] "italian" "regional"
## [199] "regional" "regional"
## [201] "regional" "regional"
## [203] "regional" "regional"
## [205] "regional" "regional"
## [207] "regional" "regional"
## [209] "regional" NA
## [211] "regional" "regional"
## [213] "regional" "regional"
## [215] "regional" "regional"
## [217] "regional" "regional"
## [219] "regional" "regional"
## [221] "regional" "regional"
## [223] "regional" "regional"
## [225] "regional" "regional"
## [227] "regional" "regional"
## [229] "regional" "regional"
## [231] "regional" "regional"
## [233] "regional" "kebab;grill"
## [235] "indian" NA
## [237] "indian,_japanese" "italian"
## [239] "italian_pizza;italian" NA
## [241] NA "regional"
## [243] NA NA
## [245] NA NA
## [247] NA NA
## [249] NA NA
## [251] NA NA
## [253] NA NA
## [255] NA NA
## [257] NA NA
## [259] NA NA
## [261] NA NA
## [263] NA "indian"
## [265] NA NA
## [267] NA NA
## [269] NA NA
## [271] NA NA
## [273] NA NA
## [275] NA NA
## [277] NA NA
## [279] NA NA
## [281] NA NA
## [283] NA NA
## [285] NA NA
## [287] "asian" NA
## [289] NA NA
## [291] "indian;north_indian;regional" NA
## [293] "chinese;chicken;asian;indian" NA
## [295] "indian" NA
## [297] "indian" NA
## [299] NA NA
So let us plot the restaurants as POINTs using the dat-R data we have downloaded. The cuisine attribute looks interesting; let us colour the POINT based on the cuisine offered at that restaurant.
dat_R <- dat_R %>%
drop_na(cuisine) %>% # Knock off nondescript restaurants
# Some have more than on classification ;-()
# Separated by semicolon or comma, so....
separate(col = cuisine, into = c("cuisine", NA, NA), sep = ";") %>%
separate(col = cuisine, into = c("cuisine", NA, NA), sep = ",")
## Warning: Expected 3 pieces. Additional pieces discarded in 1 rows [150].
## Warning: Expected 3 pieces. Missing pieces filled with `NA` in 150 rows [1, 2,
## 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, ...].
## Warning: Expected 3 pieces. Missing pieces filled with `NA` in 152 rows [1, 2,
## 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, ...].
# Finally good food?
dat_R$cuisine
## [1] "indian" "indian" "chinese" "regional"
## [5] "Mughlai" "chinese" "north indian" "chinese"
## [9] "chinese" "indian" "tandoor" "regional"
## [13] "Naga" "chinese" "chinese" "indian"
## [17] "spanish" "indian" "chinese" "indian"
## [21] "indian" "indian" "ice_cream" "chinese"
## [25] "indian" "ice_cream" "indian" "indian"
## [29] "indian" "international" "international" "ice_cream"
## [33] "seafood" "indian" "international" "chinese"
## [37] "indian" "indian" "chinese" "indian"
## [41] "international" "ice_cream" "ice_cream" "italian"
## [45] "indian" "chinese" "ice_cream" "international"
## [49] "indian" "coffee_shop" "international" "indian"
## [53] "international" "international" "international" "international"
## [57] "indian" "indian" "international" "chinese"
## [61] "indian" "indian" "chinese" "indian"
## [65] "indian" "indian" "international" "indian"
## [69] "italian" "indian" "indian" "regional"
## [73] "indian" "international" "regional" "chinese"
## [77] "pizza" "pizza" "Multi-cuisne" "indian"
## [81] "indian" "indian" "andhra" "regional"
## [85] "regional" "regional" "regional" "regional"
## [89] "regional" "regional" "regional" "regional"
## [93] "regional" "regional" "regional" "regional"
## [97] "regional" "regional" "regional" "regional"
## [101] "regional" "regional" "italian" "regional"
## [105] "italian" "regional" "regional" "regional"
## [109] "regional" "regional" "regional" "regional"
## [113] "regional" "regional" "regional" "regional"
## [117] "regional" "regional" "regional" "regional"
## [121] "regional" "regional" "regional" "regional"
## [125] "regional" "regional" "regional" "regional"
## [129] "regional" "regional" "regional" "regional"
## [133] "regional" "regional" "regional" "regional"
## [137] "regional" "regional" "regional" "regional"
## [141] "kebab" "indian" "indian" "italian"
## [145] "italian_pizza" "regional" "indian" "asian"
## [149] "indian" "chinese" "indian" "indian"
Now let’s plot the Restaurants as POINTs:
# http://www.stat.columbia.edu/~tzheng/files/Rcolor.pdf
#
ggplot() +
geom_sf(data = dat_B, colour = "burlywood1") +
geom_sf(data = dat_H, colour = "gray80") +
geom_sf(data = dat_R %>% drop_na(cuisine), aes(fill = cuisine), colour = "black", shape = 21) +
theme(legend.position = "right") +
labs(title = "Restaurants in South Central Bangalore",
caption = "Based on osmdata")
We could have done a (much!) better job, by combining cuisines into simpler and fewer categories, but that is for another day!!
Let’s go meet the Gryphon and the Mock Turtle…carrying this Pig!
rnaturalearth and tmapdata("World")
data("metro")
head(metro, n = 3)
## Simple feature collection with 3 features and 12 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 3.04197 ymin: -8.83682 xmax: 69.17246 ymax: 36.7525
## Geodetic CRS: WGS 84
## name name_long iso_a3 pop1950 pop1960 pop1970 pop1980 pop1990
## 2 Kabul Kabul AFG 170784 285352 471891 977824 1549320
## 8 Algiers El Djazair (Algiers) DZA 516450 871636 1281127 1621442 1797068
## 13 Luanda Luanda AGO 138413 219427 459225 771349 1390240
## pop2000 pop2010 pop2020 pop2030 geometry
## 2 2401109 3722320 5721697 8279607 POINT (69.17246 34.52889)
## 8 2140577 2432023 2835218 3404575 POINT (3.04197 36.7525)
## 13 2591388 4508434 6836849 10428756 POINT (13.23432 -8.83682)
tmap_mode("plot")
## tmap mode set to plotting
# Group 1
tm_shape(World) +
tm_polygons("HPI") +
# Group 2
tm_shape(metro) +
tm_bubbles(size = "pop2020", col = "red")
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(World) +
tm_polygons("HPI") +
tm_shape(metro) +
tm_bubbles(size = "pop2020", col = "red") +
tm_tiles("Stamen.TonerLabels")
## Legend for symbol sizes not available in view mode.
india <-
ne_states(country = "india",
geounit = "india",
returnclass = "sf")
india_neighbours <-
ne_states(country = (c("india", "sri lanka", "pakistan",
"afghanistan", "nepal","bangladesh")
)
)
names(india)
## [1] "featurecla" "scalerank" "adm1_code" "diss_me" "iso_3166_2"
## [6] "wikipedia" "iso_a2" "adm0_sr" "name" "name_alt"
## [11] "name_local" "type" "type_en" "code_local" "code_hasc"
## [16] "note" "hasc_maybe" "region" "region_cod" "provnum_ne"
## [21] "gadm_level" "check_me" "datarank" "abbrev" "postal"
## [26] "area_sqkm" "sameascity" "labelrank" "name_len" "mapcolor9"
## [31] "mapcolor13" "fips" "fips_alt" "woe_id" "woe_label"
## [36] "woe_name" "latitude" "longitude" "sov_a3" "adm0_a3"
## [41] "adm0_label" "admin" "geonunit" "gu_a3" "gn_id"
## [46] "gn_name" "gns_id" "gns_name" "gn_level" "gn_region"
## [51] "gn_a1_code" "region_sub" "sub_code" "gns_level" "gns_lang"
## [56] "gns_adm1" "gns_region" "min_label" "max_label" "min_zoom"
## [61] "wikidataid" "name_ar" "name_bn" "name_de" "name_en"
## [66] "name_es" "name_fr" "name_el" "name_hi" "name_hu"
## [71] "name_id" "name_it" "name_ja" "name_ko" "name_nl"
## [76] "name_pl" "name_pt" "name_ru" "name_sv" "name_tr"
## [81] "name_vi" "name_zh" "ne_id" "geometry"
names(india_neighbours)
## [1] "featurecla" "scalerank" "adm1_code" "diss_me" "iso_3166_2"
## [6] "wikipedia" "iso_a2" "adm0_sr" "name" "name_alt"
## [11] "name_local" "type" "type_en" "code_local" "code_hasc"
## [16] "note" "hasc_maybe" "region" "region_cod" "provnum_ne"
## [21] "gadm_level" "check_me" "datarank" "abbrev" "postal"
## [26] "area_sqkm" "sameascity" "labelrank" "name_len" "mapcolor9"
## [31] "mapcolor13" "fips" "fips_alt" "woe_id" "woe_label"
## [36] "woe_name" "latitude" "longitude" "sov_a3" "adm0_a3"
## [41] "adm0_label" "admin" "geonunit" "gu_a3" "gn_id"
## [46] "gn_name" "gns_id" "gns_name" "gn_level" "gn_region"
## [51] "gn_a1_code" "region_sub" "sub_code" "gns_level" "gns_lang"
## [56] "gns_adm1" "gns_region" "min_label" "max_label" "min_zoom"
## [61] "wikidataid" "name_ar" "name_bn" "name_de" "name_en"
## [66] "name_es" "name_fr" "name_el" "name_hi" "name_hu"
## [71] "name_id" "name_it" "name_ja" "name_ko" "name_nl"
## [76] "name_pl" "name_pt" "name_ru" "name_sv" "name_tr"
## [81] "name_vi" "name_zh" "ne_id"
#ind <- metro$iso_a3 == "IND"
#metro_ind1 <- metro[ind,]
metro_ind2 <- subset(metro, iso_a3 == "IND")
metro_neighbours <- metro %>% dplyr::filter(iso_a3 %in% c("IND","PAK", "LKA", "BGD","NPL"))
tm_shape(World %>% dplyr::filter(iso_a3 %in% c("IND", "AFG", "PAK", "NPL", "BGD", "LKA"))) +
tm_borders() +
tm_shape(india) +
tm_polygons("name") +
tm_shape(metro_ind2)+
tm_dots(size = "pop2020") +
tm_layout(legend.outside = TRUE,legend.outside.position = "right") +
tm_credits("Geographical Boundaries are not accurate",size = 0.5,position = "right") +
tm_compass(position = c("right", "top")) +
tm_scale_bar(position = "left") +
tmap_style("watercolor") #cobalt #gray #white #col_blind #beaver #classic #watercolor #albatross #bw
## tmap style set to "watercolor"
## other available styles are: "white", "gray", "natural", "cobalt", "col_blind", "albatross", "beaver", "bw", "classic"
## Credits not supported in view mode.
## Compass not supported in view mode.
## Warning: Number of levels of the variable "name" is 35, which is
## larger than max.categories (which is 30), so levels are combined. Set
## tmap_options(max.categories = 35) in the layer function to show all levels.
## Legend for symbol sizes not available in view mode.
tm_shape(india_neighbours) +
tm_polygons("name") +
tm_shape(metro_neighbours) +
tm_dots(size = "pop2020") +
tm_layout(legend.outside = TRUE,legend.outside.position = "right") +
tmap_options(max.categories = 10) +
tm_credits("Geographical Boundaries are not accurate",size = 0.5,position = "center")
## Credits not supported in view mode.
## Warning: Number of levels of the variable "name" is 112, which is
## larger than max.categories (which is 10), so levels are combined. Set
## tmap_options(max.categories = 112) in the layer function to show all levels.
## Legend for symbol sizes not available in view mode.
tmap_mode("plot")
## tmap mode set to plotting
tm_basemap("Stamen.Watercolor") +
tm_shape(metro, bbox = "India") +
tm_dots(col = "red",
# user-chosen group name for layers
group = "Metropolitan Areas") +
tm_shape(World) +
tm_borders() +
tm_tiles(server = "Stamen.TonerLabels",group = "Labels") + # ADDS LABELS!!!
tm_graticules()
Emine Fidan, Guide to Creating Interactive Maps in R
Nikita Voevodin,R, Not the Best Practices
Draw a map of your home-town with your favourite restaurants shown. Pop-ups for each restaurant will win bonus points.